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Abstract. Although QR iterations dominate in eigenvalue computations, there are several im- 
portant cases when alternative LR-type algorithms may be preferable. In particular, in the symmetric 
tridiagonal case where differential qd algorithm with shifts (dqds) proposed by Fernando and Parlett 
enjoys often faster convergence while preserving high relative accuracy. In eigenvalue computations 
for rank-structured matrices QR algorithm is also a popular choice since, in the symmetric case, 
the rank structure is preserved. In the unsymmetric case, however, QR algorithm destroys the rank 
{^JQf structure and, hence, LR-type algorithms come to play once again. In the current paper we adapt 

i—{ ' several variants of qd algorithms to quasiseparable matrices. Remarkably, one of them, when applied 

to Hessenberg matrices, becomes a direct generalization of dqds algorithm for tridiagonal matrices. 
Therefore, it can be applied to such important matrices as companion and confederate, and pro- 
vides an alternative algorithm for finding roots of a polynomial represented in a basis of orthogonal 
polynomials. Results of preliminary numerical experiments are presented. 
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1. Introduction. Eigenvalue problem for rank-structured matrices such as se- 
miseparable, quasiseparable, unitary-plus-rank-one and others has been an area of 
intense research in the last decade. This is due to the fact that the class of rank- 
structured matrices includes many other important classes, among which are banded 
matrices and their inverses, unitary Hessenberg matrices, companion and confederate 
matrices. Moreover, computations with such matrices are cheap and one step of any 
iterative algorithm usually takes only 0{n) arithmetic operations, where n is the size 
of the matrix. 

There is one significant difference of rank-structured eigenvalue computations 
from the unstructured case. A prospective algorithm must preserve the low-rank 
structure of the initial matrix (maybe in some other form) to take advantage of fast 
linear algebra. This conservation property may not always be taken for granted and 
must be taken into account while developing new algorithms. 

As usual, there are two competitive approaches to eigenvalue "hunting" : QR- 
and LR-type algorithms. Let us summarize the current state of the above mentioned 
approaches. 

QR: Development of QR-type eigenvalue solvers for rank-structured matrices was 
initially motivated by the matrix formulation of polynomial root-finding prob- 
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This approach was successfuUy pursued in [51[7] where certain low-rank preser- 
vation properties of QR-itcrations for companion matrices were used. It was 
observed soon [BJ |51 13] that more gencrahy QR iterations preserve the struc- 
ture of unitary-plus-rank-one matrices and companion matrix is only one 
particular representative of this class. Another direction of research aims at 
finding eigenvalues of symmetric quasiseparable matrices |151 123j , the struc- 
ture that is also preserved under QR iterations. 
LR: One clear advantage of LR-type algorithms over the QR-type ones is that the 
former preserve quasiseparable structure of iterates even in the unsymmetric 
case. The disadvantage is the possible instability in floating point arithmetic. 
However, as it has been observed recently, even the use of orthogonal trans- 
formations for rank-structured matrices may lead to an unstable algorithm 
[T^ . There are two papers to be mentioned here: [H] gives a Cholesky LR 
algorithm for the symmetric positive definite quasiseparable matrices and |3] 
uses Neville representation of quasiseparable matrices to develop a qd-type 
method. 

In the important case of symmetric positive definite tridiagonal matrices Fernando 
and Parlett [161 US developed a root- free eigenvalue problem solver called differential 
qd algorithm with shifts, or dqds for short. The algorithm is fast and accurate and 
has become one of LAPACK's main eigenvalue routines. Although [3] attempted to 
transfer this algorithm to the quasiseparable case, the proposed algorithm is not dqds 
in the strict sense as it uses Neville-type eliminations. 

The main contribution of the current paper is a direct generalization of dqds 
algorithm of Fernando and Parlett to the quasiseparable matrix case. By achieving 
this goal we strictly follow the methodology of [20] and, as a by-product, derive 
several new eigenvalue algorithms applicable in different cases. We list below all the 
algorithms obtained in the paper. 

• New LU decomposition algorithm for general quasiseparable matrices. 

• Stationary and Progressive qd algorithms with shifts for general quasisepa- 
rable matrices. 

• Differential qd algorithm with shifts (dqds) for Hessenbcrg quasiseparable 
matrices. 

• Specification of dqds algorithm for companion and comrade matrices. 

In the symmetric positive-definite tridiagonal case both dqds algorithm of Fer- 
nando and Parlett and the variant of QR algorithm by Demmel and Kahan [llj 
guarantee high relative accuracy of the computed eigenvalues (although dqds algo- 
rithm is not backward stable). The unsymmetric case is much more complex and, 
to the best of our knowledge, relative accuracy has not been proved for any of the 
existing algorithms even in the simple tridiagonal case. In this paper we are mainly 
interested in the unsymmetric eigenvalue problem and do not address the issue of rel- 
ative accuracy. However, we can say definitely that the proposed dqds algorithm for 
Hessenberg quasiseparable matrices is not backward stable as it is a straightforward 
generalization of tridiagonal dqds algorithm. Another issue is that the new algorithm 
uses LU factorization without pivoting at the initial step. We, therefore, assume that 
this factorization exists. Nevertheless, numerical experiments with many different 
matrices show that the new dqds algorithm often delivers more accurate result than 
its QR-bascd competitors. 

The outline of the paper is as follows. In Section [2| we recall the definition of 
quasiseparable matrices and also derive their LU factorization that will be our main 
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tool later. A restricted version of this algorithm applicable to diagonal plus semise- 
parable matrices (a subclass of quasiseparable matrices) was derived in [T7] and some 
formulae of it are similar to the ones in the inversion algorithm given in |13j . Still, new 
LU factorization algorithm, to the best of our knowledge, has never been published 
and may be useful to those who need a fast system solver for quasiseparable matrices. 
The new algorithm uses the idea of successive Schur complements and in this respect 
it is also different from the algorithm proposed in [22l p. 171] that is based on the 
representation of the original matrix as a sum of lower and upper triangular matrices. 
Moreover, we provide an explicit pseudocode of the algorithm, while no pseudocode 
was provided in [22]. The complexity of the algorithm is 0{N) and it is valid in the 
strongly regular case (i.e. all its block leading principal minors are non-vanishing) . In 
the subsequent Section [3] we present two versions of qd algorithm for general quasi- 
separable matrices: stationary and progressive. Section 2] is central in the paper and 
the new dqds algorithm is presented there. Section [5] provides an alternative deriva- 
tion of dqds algorithm via the generalized Gram-Schmidt process and also shed light 
on the meaning of some parameters arising in the algorithm. We next specialize the 
dqds algorithm for general Hessenberg quasiseparable matrices to more specific cases 
of companion and comrade matrices in Section |6l Results of preliminary numerical 
experiments are presented in the final Section [T] 

2. Quasiseparable matrices. There arc many important classes of structured 
matrices with the property of having low-rank blocks above and below the diagonal 
that one can meet in applications. Among the most well-known are semiseparable, 
quasiseparable, 77-matrices, TJ^-matrices and mosaic-skeleton matrices. In the current 
paper we are particularly interested in eigenvalue problem for quasiseparable matrices 
leaving possible extensions of the results to other classes of rank-structured matrices 
for future research. We will also consider the most general version of quasiseparable 
matrices, namely block quasiseparable matrices as many of our results trivially gener- 
alize from the scalar to the block case. The usual definition of a block quasiseparable 
matrix is given below. 

Definition 2.1 (Rank definition of a block quasiseparable matrix). Let A he a 
block matrix of block sizes {nk}^^i then it is called block {r\r^)- quasiseparable if 



where r' (r^) is called lower (upper) order of quasiseparability. 

In other words, quasiseparable matrices are those having low rank partitions in 
their upper and lower parts. In what follows we will call block {r\r")- quasiseparable 
matrices simply quasiseparable for shortness. 

In order to exploit the quasiseparability of matrices in practice one must use a 
low-parametric representation of them. There are many alternative parametrizations 
of quasiseparable matrices all of which use 0{N) parameters, where N is the size of the 
matrix. Having such a parametrization at hand one can write most of the algorithms, 
e.g., inversion, LU, QR, matrix- vector multiplication in terms of these parameters 
and, therefore, the complexity of these algorithms is 0{N). In the current paper 
we will use the so-called generator representation (Definition 12.21 below) proposed by 
Eidelman and Gohberg [T3] . For alternative parametrizations see [23 (TU] . 

Definition 2.2 (Generator definition of a block quasiseparable matrix). Let A 



max rank A{K + 1 : N,l : K) < r' 



max rank A(l : K,K + 1 : N) < r". 
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be a block matrix of block sizes {nkYj^^i then it is called block {r\r'^)-quasiseparable 
if it can be represented in the form 



di 

P2qi 

P3a2qi 



gih2 
d2 

Pzq2 



.91^2/13 
52/13 
dz 



Pnan-l ■ ■ ■ a2qi Pna„-1 . . . 0392 PnOn-l 



5162 • • • bn-lKi 
52^3 • • ■ K-lhn 

gsbi . . . bn-ihn 



(2.1) 



where parameters (called generators) {dk, qk, Oik, Pk, gk, bk, hk} are matrices of sizes 
as in the table below. 

Table 2.1 

Sizes of generators of a block quasiseparable matrix. 



dk 



qk 



a-k 



Pk 



Qk 



# of rows Uk 



nk 



nk 



'k-i 



'k^i 



# of cols Uk Uk 



' k-i 



'k-i 



r'k 



Orders of quasiseparability r and r " are maxima over the corresponding sizes of 
generators: 

r' = max ri, r" = max rlf. 

l<fc<n-l l<fc<n-l 



Remark 2.3. Generators are not unique, there are infinitely many ways to 
represent the same quasiseparable matrix with different generators. For the relation 
between different sets of generators see \14^ . 

Remark 2.4. Sizes r\, and r'^ of generators are directly related to ranks of sub- 
matrices in the lower and upper parts correspondingly. Namely, if we let K to be the 
k 

block index: A" = Ui, then 
1=1 

T:m\kA{K + l:N,l:K)<r{, rank A(l : A', AT + 1 : iV) < r^. (2.2) 

Moreover, for any quasiseparable matrix there exists a set of generators with sizes r\, 
and r^ that satisfy (|2.2p with exact equalities (such generators are called minimal^. 
For their existence and construction see m^ . 

We next derive LU factorization algorithm for a general block quasiseparable 
matrix in terms of the generators it is described by. First, let us note that quasise- 
parable structure of the original matrix implies the quasiseparable structure of L and 
U factors. The theorem below makes this statement precise and, in addition, relates 
generators of an original matrix to the generators of its factors. 

Theorem 2.5. Let A be a strongly regular N x N block {r^r"^)- quasiseparable 
matrix given by generators {dk,qk,ak,Pk,gk,bk,hk} as in Let A = LU be its 

block L U decomposition of the same block sizes. Then 

(i) Factors L and U arc (r',0)- and {0,r^)-quasiseparable. Moreover, r'j,{L) = 

rl{A) and r]^{U) = r{{A) for all k = 1, . . . ,n - 1. 
(ii) L and U are parametrized by the generators {Ik, qk,ak,Pk, 0,0,0} and {dk, 
0, 0, 0, gk, bk, hk}, where Ik are identity matrices of sizes Uk x Uk and new 
parameters qk, dk and gk can be computed using the following algorithm: 
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Algorithm 1 Fast quasiseparable LU decomposition. 
Input: dk,qk,ak,Pk,^k,bk,hk 

di^di, q^=qid^^, 9i ^ 91, /i = 

foi^fc = 2 to n 1 do 

dk = dk - Pkfk-ihk 
Qk = {qk - akfk-ihk)dl^ 
gk = 9k -Pkfk-ih 
fk = akfk-ibk + qkHk 
end for 

dn dyi Pjifn—lhn- 

Output: dk,qk,gk 



Proof. Statement (i) of the theorem follows from statement (ii), so we only need 
to prove the latter part. 

fe 

Denote, as usual, by K the block index: K ~ rii and note that quasiseparable 

i=i 

representation (j2.ip implies the following recurrence relation between the blocks of A: 



A = 


■ A{1 : K, 


l:K) 


GkHk+i 


Pk+l 


Qk 


A{K + 


I: N,K ^ 


~ 1 : 


N) J 


Qi - 


qi, Qk -- 


= [akQk 


-1 qk]. 


k = 2,.. 


n - 


-1; 


Pn = 


Pn, Pk = 


= [Pk ; Pk+iak], 


k = n — 


1,. 


..2; 


Gi = 


91, Gk = 


= [Gk-i 


bk ; 9k], 


k = 2,. 


. n 


-1; 


Hn = 


hn , Hk 


= [hk bkHk+i], 


k = n ~ 


1,. 


..2. 



(2.3) 



The proof will be constructed by induction. We will show that for each k 



Ak 
^11 


GkHk+i 




^11 


■ 




Tjk 


GkHk+1 


Pk+lQk 






Pk+lQk 









-k 


= 1 we get from (|2.4|): 
















di = di, 






di = 


= a 


1, 






P2Q1 = PiQid 


1, 


qi -- 


= 910^1 \ 





(2.4) 



G1H2 = G1H2, 



91 = 91- 



Let us introduce an auxiliary parameter fk = QkGk- It is easy to show by using 
2.3p that fk satisfies the recurrence relation 

fi = qi9i, fk = akfk-ibk + gfcfffc- (2.5) 
Assume that ()2.4p holds for some fixed fc, then it holds for fc + 1 if 

dk+i = [Pk+iQk 1] • [Gfc/ife+i ; dk+i], (2.6) 

Pk+2qk+l = Pk+2Qk+l ' [Gfc/lfe+l ; dk+l], (2.7) 

9k+iHk+2 = [pk+iQk 1] • Gk+iHk+2. (2.8) 
The first equality simplifies 

dk+i " Pk+iQkGkhk+i + dk+i = Pk+ifkhk+i + dk+i- 
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The second equality (j2.7p holds if 



qk+i = Qk+i ■ [Gkhk+1 ; dk+i] = 

= [a/c+iQfe Ik+i] ■ [Gkhk+i ; dfc+i] = a-k+ifkhk+i + qk+idk+i- 
Finally, the last equality (|2.8p is true if 

9k+i = [vk+iQk 1] • Gk+1 = [pk+iQk 1] • [Gkbk+1 ; fffc+i] = Pk+ifkbk+i + gk+i- 

Matrix dk+i is invertible because, by our assumption, matrix A is strongly regular. 
Hence, we conclude that (|2.4p is true also for index fc + 1 if generators qk+i, dfc+i and 
Qk+i are those computed in Algorithm [TJ The assertion of the theorem holds by 
induction. 



3. qd algorithms for general quasiseparable matrices. In the previous 
section we have shown that factors in the quasiseparable LU factorization retain the 
low-rank structure of the original matrix. It turns out that even stronger result can 
be proved, namely that the low-rank structure is preserved under iterations of the LR 
algorithm. For completeness let us prove this simple theorem below. 

Theorem 3.1. Let A be anNxN block {r^r"")- quasiseparable matrix and A — al 
be strongly regular. Define one step of shifted LR iterations by 

A-al = LU, A' = UL. 

Then A' is a strongly regular block {r'' jr"^)- quasiseparable matrix. 

Proof. From Theorem 12. 51 we know that matrices L and U are (r',0)- and (0,r")- 

k 

quasiseparable. Let {nk}]t^^ be sizes of blocks of A and let = ^ n.j. Then for 

1=1 

each K let us write the product of U and L terms in the block form assuming that 
A'li = A'{1 :K,1: K): 



A' 
^11 


A' ' 

^12 




■ C/ii 


C/l2 " 


A' 

^21 


A' 

^22 . 







U22 \ 








L21 


-^22 



Hence, 



rank(A2i) = rank(L2iC^ii) = rank(L2i) = r\ 
rank(A'^2) = i'ank(J7i2L22) = rank([/i2) 



where Un is invertible by our assumption on the strong regularity of ^ — aL . □ 

The assertion of Thcorcm l3. li lies in the heart of fast LR-typc algorithms proposed 
in [3T] and [3] . In fact Algorithm [1] can be used to derive a new version of quasise- 
parable LR algorithm but we will not do it here as our main interest lies in deriving 
qd-type algorithms that are believed to be better in practice. Below we derive 3 new 
algorithms: stationary qd, progressive qd and Hessenberg dqds. Our ultimate goal is 
the last algorithm and the first two can be seen as intermediate results. 

3.1. Stationary qd algorithm. Triangular factors change in a complicated way 
under translation. Let L and U be quasiseparable factors as in Theorem [231 our task 
is to compute L and U so that 



A - al = LU - al = LU 
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for a given shift a. Knowing that factors L and U on the left as well as L and U on 
the right are quasiseparable matrices we want to find formulae that define the direct 
mapping from the generators of the first pair to the generators of the second pair. Let 
factors L and U be given by generators {g^, Ofc, pk \ and {c?fe, g^, 6^, hk\, respectively. 
Simply inverting^ the formulae in Algorithm [1] to get quasiseparable generators of 
LU — (tI and LU and equating them we obtain: 

flfe = Ofe, bk = bk, Pk = Pk, hk = hk, 
fi = 0, fk+i = akfkbk + qk9k, fi = 0, fk+i = akfkbk + qk9k, 
akfkhk + Qkdk = akfkhk + qkdk, Pkfkbk + Qk = Pkfkbk + 
Pkfkhk + dk- cflk = PkJkhk + dk- 

The above written formulae give rise to the algorithm of computing generators {git, 
Ofc, Pfc} and {dk, gk, bk, hk} of the new factors L and U, respectively. 



Algorithm 2 Stationary qd algorithm with shift (stqd(CT)). 

Ingut: dk, qk, ak, Pk, Qk, bk,Jik and a 
di = di - all, qi = qidid^^, gi = gi 
for fc = 2 to n 1 do 

tk = ak-itk-ibk-i + qk-Wk-i - qk-idk-i 
dk = Pktkhk + dk - alk 

qk = {aktkhk + qkdk) \ gk = Pktkbk + gk 

a-k = ak, Pk = Pk, bk = bk, hk = hk 
end for 

tn = 0,n-ltn-lbn-l + 9,1-1.971-1 — 9n-1.9n-l 
dn Pn^nhn ^" dn a In, ^Pn^ Pn, 

Output: dk, qk, ak, Pk, dk, bk, hk 



Let us show that Algorithm [2] is the direct generalization of stationary qd algo- 
rithm for tridiagonal matrices. L and U factors of a normalizecQ tridiagonal matrix 
with subdiagonal entries Ik and diagonal entries Uk have very special sets of quasise- 
parable generators, namely {1, Ik, 0, 1, 0, 0, 0} and {uk, 0, 0, 0, 1, 0, 1}. For these 
special generators formulae in Algorithm [5] reduce to 

Uk = dk^ Pktkhk + dk - a ^ qk-i + dk - (t - qk-i = h-i + Uk - a - h-i, 
Ik = Pk+iqk = Pk+i {aktkhk + qkdk) d^^ = {pk+iqk)dkd^^ = IkUk/uk- 

Equation p.ip for Ik and Uk in terms of Ik and Uk is exactly the one defining stationary 
qd algorithm with shift for tridiagonal matrices (see [20l page 465]). 

3.2. Progressive qd algorithm, hi deriving progressive qd algorithm with 
shift we seek for the triangular factorization of A' — al: 

A' - cri = UL- crl = LU (3.2) 

for a suitable shift a. Theorem 13.11 tells us that UL — <tI is a quasiseparable matrix 
of the same order as L and U factors and that suitable quasiseparable generators 



all entries in positions + I) are ones 
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paramctrizatioii of L and U exists. However, we cannot obtain generators of UL — 
cri by simply inverting formulae in Algorithm [T] as was done in the derivation of 
stationary qd algorithm. For this purpose we need to exploit the quasiseparable 
matrices multiplication Algorithm 4.3 from [14] that defines generators of the product 
of two quasiseparable matrices. 

Let {Ik, qk, ak, Pk, 0, 0, 0} and {dk, 0, 0, 0, gk, bk, hk} be generators of factors 
L and U in (|3.2p . then generators {dj., q'f., aj., p'^, g^, 5'^,, /ij.} of A' — al = UL — al 
are as follows: 

<lk = 1k, g'k'^gk, a'k^ak, b'f.^bk, 

fn = 0, fk-l = bkfkClk + hkPk, 

p'k ^ dkPk + gkfkCik, K^hk + bkfkQk, 
d'k ^ dk+ gkfkQk - cr/fc. 

Equating generators in (|3.3p to the generators of LU we obtain progressive qd 
algorithm with shift described below. 



Algorithm 3 Progressive qd algorithm with shift (qds((T)). 
Input: rffe, qk^ ak, Pk, gk, bk, hk and a 

Pn — Pn , hji = hn , fn ~ 

for fc = n — 1 to 2 do 

fk ~ bk+lfk+l^k+l + hk+lPk+1 

Pk = dkPk + gkfkO-k, hk = hk + bkfkqk 
ak = ak, bk = bk 
end for 

/i = b2f2a2 + /12P2, di = di + g\f\q\ - oh 
qi = qid^^, gi=gi, /i = 
for fc = 2 to n — 1 do 

fk = a.k-lfk-lbk-1 + qk-igk-l 

qk = {qk - akfkhk)dl^, gk = gk - Pkfkbk 

dk = dk + gkfkqk - <ylk - Pkfkhk 
end for 

fn = an-lfn-lbn-1 + <Zn-l,9ri-l7 dn = dn — (jIn — Pnfnhn 

Output: dk, qk, ak, Pk, gk, bk, hk 



In the case of normalized tridiagonal matrices, L and U factors have generators 
{1, Ik, 0, 1, 0, 0, 0} and {uk, 0, 0, 0, 1, 0, 1}, and Algorithm [3] simplifies to 

h ^ Pk+iqk = {dk+iPk+i){qk/dk) ^ IkUk+i/uk, 

Uk = dk = dk +Pk+iqk - cr - Pkqk-i = Uk + Ik - u - Ik-i, 

where by Ik /Ik and Uk/uk we denoted subdiagonal entries oi L/ L factor and diagonal 
entries oi U/U factor correspondingly. One can easily spot in p.4p the progressive qd 
algorithm with shift for tridiagonal matrices [20l page 467]. 

4. Differential qd algorithm in the Hessenberg quasiseparable case. Fer- 
nando and Parlett |16| I20j proposed a modified version of progressive qd algorithm 
called differential qd algorithm. The proposed algorithm has received major atten- 
tion due to its accuracy and speed, and is implemented as DLASQ in LAPACK. It 
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turns out that differential version of progressive qd algorithm does not exist in the 
general case (the reason for this will be given in Section [S]). However, we found that 
it exists for Hessenberg quasiseparable matrices that is an important case on its own 
as matrices related to polynomials, e.g. companion, confederate, fellow, are usually 
Hessenberg. 

Consider a general Hessenberg quasiseparable matrix: 







gih2 5i&2/i3 
S2 ds 



gib2 ■ ..hn-iK 

3263 • ■ • bn-lK 

gsbi . . . bn-ih„ 



(4.1) 



dn 



Because of the Hessenberg structure, quasiseparable generators of the lower part are 
special, namely qk = Sfc, Cfc = and pk = 1. Let us note that from now on we restrict 
matrices to the scalar case, i.e. generators such as dk are all scalars. 

We next obtain differential version of Algorithm [S] Our initial derivation may 
seem cumbersome and very technical but a much cleaner and justified derivation will 
be given later in Section [S] 

First, observe that due to the special Hessenberg generators many formulae in 
Algorithm |3] can be simplified. In particular fk = hk+i and fk = Qk-idk-i- We next 
change the way dk are computed: 

dk = dk + SkSkhk+i - Sk-iSk-ihk - cf 

= dk + Sk{gk + Sk-igk-ibk)hk+i - Sk-ifjk-ihk - o 

= dk + Skdkhk+i - Sk-igk-i(hk - Skbkhk+i) - a 
= dk + Skgkhk+i - Sk-igk-ihk - cr. 

Next, let us define an auxiliary variable: 

tk = dk - Sk-idk-ihk - (J dk - Skdkhk+i)- (4.2) 

Observe that 

tk+i = dk+i - {skdk+i/dk)gkhk+i - cr = 

= dk+i/dk{dk - Skdkhk+i) - cr = dk+itk/dk - a. (4.3) 



Using identities (|4.2p and (|4.3p we can now derive dqds algorithm for Hessenberg 
quasiseparable matrices. 
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Algorithm 4 Differential qd algorithm with shift (dqds(cr)), Hesscnberg case. 
Input: Sfc, dfc, gk, 6fe, hk and a 

ti=di-(T^ ?i=.9i, di = + sigi/i2, 

si^ sid2/di t2 =tid2/di - a 

for fc = 2 to n — 1 do 
hk = hk + Skbkhk+i 

dk = gk - Sk-idk-ibk 

dk = tk + Skdkhk+i 
Sk = Skdk+i/dk^ 
tk+i = tkdk+i/dk - a 
end for 

dn tn 7 

Output: Sk, dk, gk, hk 



Let us show that Algorithm 2] trivially reduces to the well-known dqds algorithm 
in the normalized tridiagonal case. In this case we have bk = and gkhk+i = 1 for 
all k and, hence, 

dk = tk + Sk, 

Sk = Skdk+i/dk, 

tk+i ~ tk{dk+i/dk) - cr. 

These formulae constitute dqds(cr) algorithm for normalized tridiagonal matrices, see 
1201 page 468]. 

In the tridiagonal case dqds algorithm can be run from bottom to top of L and U 
factors. Let us show how Algorithm |4] in the Hessenberg case can be also transformed 
to work backward. Let A = LU be Hessenberg quasiseparable matrix with generators 
{sk, dk, gk, bk, hk} of L and U factors. Roots of characteristic polynomial det(a;/ — A) 
equal to the roots of det(a;/— JAJ), where J is antidiagonal matrix. Transform JAJ is 
simply backward reordering of rows and columns of A. Applying this transform to the 
LU factorization we get JAJ = {JLJ){JU J). It is trivial to show that matrices {JLJ) 
and ( JUJ) retain all the properties of L and U and their quasiseparable generators are 
{s„_fc, dn-k+i, hj^_k+i, b'^_k+i, Qn-k+i}- Hence, running AlgorithmHin a forward 
way on JLJ and JU J is same as running it backward on L and U. 

5. Generalized Gram Schmidt process and quasiseparable dqds algo- 
rithm. Parlett [20] has shown that dqd algorithm for tridiagonal matrices can be 
interpreted as the generalized Gram-Schmidt orthogonalization process applied to 
matrices L and U from A' = UL. In this section we follow the same approach and 
derive Algorithm |4] directly from the Gram-Schmidt and without referencing to Al- 
gorithms |3l Our analysis uses little of quasiseparable matrices theory and is very 
accessible. 

Let F = [/i /2 ■ • ■ fk] and G = [gi g2 ■ ■ ■ gk] be a pair of vector sets. 
Next theorem relates biorthogonalization of F and G (generalized Gram-Schmidt) to 
the LU factorization of G*F (here and thereafter A* denotes the conjugate transpose 
of A). 

Theorem 5.1 (Theorem 1, 2^). Let F and G he complex nxk matrices, n > k, 
such that G*F permits triangular factorization: 

G*F = LDR, 
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where L and R are unit triangular (left and right), respectively, and D is diagonal. 
Then there exist unique n x k matrices Q and P such that 

F = QR, G^PL\ P*Q = D. 



Remark 5.2. In practice, when n ^ k and D is invertible one can omit Q 
and write F = {P*)^^{DR) = {P*)~^U, G = PL* and still call it Gram-Schmidt 
factorization. The important feature is the uniqueness of Q and P. The columns of Q 
and rows of P* form a pair of dual bases for the space of n-vectors (columns) and its 
dual (the row n-vectors). There is no notion of orthogonality or inner product here; 
p*qj = simply says that p* annihilates qj, i =^ j . 

To derive the unshifted version of dqds Algorithm |4] apply Gram-Schmidt to the 
columns of L and U*, in the natural order, to obtain L and U, then according to 
Theorem O 

UL = LU. 

Since matrix L is unit lower bidiagonal, the matrix P [P* in Theorem 15. ip such 
that L = P^^U and U = LP is upper Hessenberg. Fortunately, as will be shown 
soon, matrix P has a simple 0{n) representation versus dense n{n — l)/2 and makes 
the derivation of linear complexity algorithm possible. Matrices L and U will be 
transformed to U and L correspondingly by a sequence of simple transformations. 

At the start of Gram-Schmidt factorization L and U factors are 



L 



1 

si 1 

S2 





di gih-i 


5162^3 • 


■ ■ ffl&2 ■ 


• bn-lhn 


u = 


d2 


52^3 


• • ff2&3 ■ 


■ bn-lhn 




ds ■ 


• • ff3&4 ■ 


■ bn-lhn 



Let P = Pn-i ■ ■ . P-iPi- Choose the first transformation Pi as 



Pi = 



di gih2 ^162/13 
-Si 1 

1 



5162 • ..bn-lK 





then 



Pi' = 



~gih2/di • gi&2/i3 

si/di di/di -si/di • 5162/13 



-1/di ■ gi&2 • ■ • bn-lhn 
-si/di ■ gib2 . ..bn-lhn 



where di ~ di+ sigi/i2. The careful reader may check that matrices Pi and Pj^ ^ are 
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such that 





' di 


91^2 


3162/13 •■• gib2 ■■ -bn-lhn 






1 




PiL = 






1 








S3 1 



and 



where 





1 













Sl 


t2 


52/13 


••■ g2b3 .. .bn-lhn 


C/Pf ^ = 






d3 


■■■ gsbi .. .bn-lhn 


di ~ di -\ 


- sigi 


h2, 


hk = 


hk + Skbkhk+1, 


si = sid2/di, 


92 


= 92 - 


-sigib2, t2 = did2/di 



(5.1) 



Note that (2 : 71, 2 : n) submatrices of PiL and UP^^ retained the form of the 
initial factors L and U . Therefore, we can continue the orthogonahzation process by 
using matrices Pk, shifted analogs of Pi. Formulae (|5.ip as well as the ones we will 
obtain recursively repeat formulae of Algorithm |3] in the case cr = 0. So, we have 
derived the dqd algorithm without referencing qd. The reason for the non-existence 
of the dqds algorithm in the general non-Hessenberg case is that if L factor is dense, 
then matrices Pk must also be dense and no recursive nested C(n) update of factor 
U is possible. 

One remarkable corollary of the Gram-Schmidt derivation of dqds algorithm is 
that it reveals the meaning of auxiliary quantities tk- 

Theorem 5.3. Let A be a strongly regular Hessenberg quasiseparable matrix of 
an arbitrary order. Then quantities tk generated by Algorithm [7] with zero shift are 
such that 



- - [A-%k, 



fc = 1, 



We omit the proof as it is identical to the proof of [501 Theorem 2] and only uses the 
fact that fc-th row and column of [Pk-i ■ ■ ■ Pi)L and U{Pk-i ■ ■ ■ Pi)~^ correspondingly 
are singletons. 

5.1. Incorporation of shift. In order to derive the shifted version of dqds one 
needs to apply Gram-Schmidt process to matrices U — aL^^ and L: 

UL- al ^ {U - aL-^)L [{U - cri"^)P"^] • [PL] ^ LU. 



At first glance the additional term —aL~^ appears to spoil the derivation of P. How- 
ever, since P = Pn-i ■ . ■ P2P1 and each of P^^ act only on two rows when multiplied 
from the right, it is not necessary to know all the entries of L^^ in advance but only 
the (fc + l,fc) entry immediately below the main diagonal. At the fc'th step of the 
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Gram-Schmidt process, the only change in the transformation in comparison to the 
unshifted case would be in the active 2x2 submatrix of {U — <jL^^){Pi . . . Pk-i)^^'- 



tk gkhk+i 
ask dk+i - (T 



1 

Sk 



—gkhk+i 
tk 





■ 1 







Sk 


tk+i_ 



This yields 



tk + Skdkhk+i = dk, as before, 
dk+i = Skdk, as before, 
-askgkhk+i + dk+itk — (Jtk = tk+idk- 



The last equation simplifies to tk+i 
version of dqds. 



dk+itk/dk — o", which is exactly the shifted 



6. dqds algorithm for companion and comrade matrices. In this Section 
we specialize dqds((T) algorithm for matrices related to polynomials, namely well- 
known companion and comrade matrices. In the first case the new algorithm will 
become a root-finding algorithm for polynomial written in monomial basis and in the 
second case — in orthogonal polynomials basis. 

6.1. dqds for companion matrix. Companion matrix has the form 



(6.1) 



and its most useful property for applications is that its characteristic polynomial P{x) 
is encoded in the entries: 





-rrin-i 


-mn-2 ■ ■ 


—mi 


-Too 




1 
















1 
















1 






P{x) = det(a;/ - C) = + mn-ix"-~^ H h mix + toq. 

Let a be an arbitrary shift, then LU factorization oi C — al is 
1 

Ho{<y) 



(6.2) 



L 



U 



Hoi-y) 



Hi[<y) 

H2{(J) 



HoM 



g„-2(a-) 

mi 

Ho(cr) 



Ho(>y) 



P(<y) 



(6.3) 



Where Hk{x) are celebrated Horner polynomials arising in the Horner rule of 
evaluating P{x): 



Hoix)^l, Hk{x)^x-Hk-i{x)+mn-k, fc = l,...,n, H,,{x) ^ P{x). (6.4) 
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Factorization (|6.3p exists for any a s.t. H^icr) ^ for all k. In particular it exists 
for CT = if m^, ^ 0, Vfc. Quasiseparable generators {sk, dk, gk, bk, hk} of i and U 
factors needed for the input to the dqds Algorithm 2] are given in Table 16.11 Let us 
note that recursion (j6.4l) may not be a good way of evaluating LU factorization due 
to the possibility of overflow, Algorithm [T] provides a better way to do this. 



Table 6.1 

Generators of L and U factors of a companion matrix. 

Sk dk gk bk hk 



Generators in Table 16.11 are special because bk = I for all k. It was shown in [I] 
that quasiseparable matrix is semiseparable if and only if there is a choice of generators 
such that 6fe = 1. Semiseparability means that the whole upper triangular part of 
matrix U is an upper part of some rank-one matrix. Generators bk are preserved 
under iterations of Algorithm 2] and, hence, semiseparable property is also preserved. 
We conclude that factors L and U stay bidiagonal and semiseparable correspondingly 
during the course of dqds iterations. 

6.2. dqds for comrade matrix. Let {?'fe(a;)}fc=o be a system of (monic) or- 
thogonal polynomials associated with some nonnegative measure on the real line M. 
As is well known, every such system satisfies a recurrence relation of the form 

rfe+i(x) = {x - ak) ■ rk{x) - Pk ■ rk-i[x), fc = 1,2,3, ... , 
ro(a:;)=0, ri{x) ^ {x ~ a^) ■ rii{x), 

where Uk and j3k > Q are unique (for a given measure and support) real constants. 
Let P{x) be a monic polynomial of degree n represented in the basis of {r^ (2;)}^_q : 

P{.x) = rn{x) + m„_i ■ rn-i{x) H V mi ■ ri{x) + rriQ ■ ro(x). 



Then comrade matrix for P{x) is defined as follows: 



C 



Q;„_i-m„_i /3„_i-m„„2 ■■■ -mo 

1 a„_2 /3„-2 



1 



ai Pi 

1 OiQ 



(6.5) 



For an arbitrary value of the shift a, the LU factorization of shifted comrade 
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matrix C — al is 



L = 



1 



Hi(<y) 

H2((J) 



U = 



Ho(a) 



Ho(cr) 



"777(57 



Ho(a) 

■777130 



_P(£) 



(6.6) 



where Hk{x) are generalized Horner polynomials defined via the recurrence relation 
(Clenshaw algorithm [3]): 



Hkix) = (x - a„_fe) • Hk-iix) ~ l3n-k+iHk-2{x) + m„_fe, fc = 2, . . . ,n, 
iJo(a;) = 1, -ffi(x) = x - a„_i + 7n„_i, Hn{x) ^ P{x). 



3.7) 



LU factorization (|6.6p exists if and only if Hk{(j) 7^ for all A; = 1, . . . , n. Matrix 
U is quasiseparable of order 2 and, therefore, its generators gk, bk, hk are matrices of 
appropriate sizes. All generators of L and U factors are given in Table [621 

Table 6.2 

Generators of L and U factors of a shifted comrade matrix. 



Sk 



9k 



bk 



Hk-i{<y) 




1 



k+1 TTLn-k 



Generators in Table [6?2] are input parameters to the dqds((T) Algorithm |4] that one 
can run with suitable shifts to compute roots of polynomial P{x). Due to the special 
form of generators of U and some conservation properties of Algorithm 21 matrix U 
stays the sum of a bidiagonal and semiseparable of order one at every step of dqds 
iterations. 

6.3. dqds for fellow matrix. One way to evaluate roots of a polynomial rep- 
resented in the basis of polynomials orthogonal on the unit circle (also called Szego 
polynomials) , is to form a so-called fellow matrix. Fellow matrix is also Hessenberg 
and is an analog of comrade matrix. Hovewer, in this case its upper triangular part 
is dense. It is known that fellow matrix is quasiseparable of order 2 (see, for instance, 
[2]) and, therefore, Algorithm |4] applies to it. We do not derive initial LU factorization 
of a fellow matrix here but, if needed, it can be done by analogy with companion and 
comrade cases. 

7. Numerical results. In this section we present results of numerical tests with 
the new dqds Algorithm H implemented in MATLAB (2010b, The MathWorks Inc., 
Natick, MA). Wc compare performances of the following algorithms for finding roots 
of polynomials: 
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• qs_dqds — qusiseparablc dqds Algorithmic 

• compan_qr — companion QR algorithm described in [i]. The implementation 
of this algorithm was generously provided by Paola Boito and is also accessible 
at www . unilim . f r/pages_perso/paola . boito/sof tware . html; 

• eig — MATLAB's main eigenvalue routine. When applied to a comrade 
matrix ()6.5p can be used to find roots of polynomials represented in orthogonal 
polynomials' bases. 

• roots — MATLAB's main routine for finding polynomial's roots. Same as eig 
applied to a companion matrix. 

Also in the tables of this section we will use the following notations: 

• n — polynomial's degree; 

• ni — number of iterations per root used by an algorithm; 

• e — relative accuracy of the computed roots evaluated as £ = max,; , 
where .t,; and Xi are exact and computed roots, respectively. 

It is a common strategy to perform diagonal scaling (balancing) of a matrix before 
applying an eigenvalue algorithm (see, for instance, |19|). Balancing can greatly 
improve accuracy of the computed eigenvalues. Unfortunately, algorithm compan_qr 
does not allow balancing because it represents initial companion matrix and all the 
subsequent iterates as unitary plus rank-one matrices and balancing destroys this 
structure. For completeness we present results of all the remaining algorithms with 
and without balancing. For balancing algorithm we simply used MATLAB's balance 
routine. 

Coefficients of test polynomials in monomials basis were computed in MATLAB 
from their roots using the poly command. Coefficients in orthogonal polynomials 
bases were computed from the corresponding coefficients in monomials basis (details 
of this algorithm can be found, for instance, in [5]). 

Details of our implementation of qs_dqds algorithm are as follows. The initial 
LU factorization was computed with zero shift. Value of the shift a at every con- 
secutive iteration was chosen to be the current iterate's (fc, k) entry (when comput- 
ing fc-th eigenvalue). We used a simple deflation criterion < £|^fc_fc| with 
e = 10"^^. Since, at the moment qs_dqds algorithm is only developed for single 
shifts, in our experiments we have considered polynomials with real roots. Current 
implementation and the script that reproduces results of this section are avalable at 
,www ■ maths . ed . ac. uk/~pzhlobic/sof tware . shtml , 

n 

Test 1: Wilkinson's first polynomial p{x) = H ^ 

It is known due J.H. Wilkinson that large roots of this polynomial are extremely 
sensitive to the perturbation of coefficients. Hence, computing roots of this polynomial 
is a hard test for any root-finding algorithm and we can't expect any such algorithm to 
compute these roots accurately for a polynomial of high degree. Results of numerical 
tests presented in Table [771] demonstrate that algorithm qs_dqds always delivers same 
or higher relative accuracy and the number of iterations per root is quite moderate (3.3 
on average). In the companion matrix case our implementation modifies parameters 
"in place" and, therefore, uses only An storage and lOn arithmetic operations per 
iteration. 

n 

Test 2: Reversed Wilkinson's first polynomial p{x) = J| (x — 1/i). 

1=1 

This polynomial happened to be a harder test for the algorithms. Roots computed 
by qs_dqds were from 10 to 10"* times more accurate than those computed by MAT- 
LAB's roots function. Also note that surprisingly balancing did not have any influence 
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Table 7.1 
Wilkinson's first polynomial. 



with balancing without balancing 





ns dnd 


s 




roots 


qs_dc|cl 


s 


compan_qr 


roots 


n 


£ 


ni 


e 


e 


ni 




e 


ni 






10 


1.6e-ll 


3 


4 


5.6e-10 


2.1e-ll 


3.3 


7 


le-11 


4.0 


5 


Oe-11 


11 


9.3e-ll 


3 


4 


2.4e-09 


7.5e-ll 


3.3 


1 


3e-10 


4.1 


1 


lc-09 


12 


2.5e-09 


3 


4 


2.7e-08 


2.4e-09 


3.3 


1 


6e-09 


4.1 


2 


9e-09 


13 


1.2e-08 


3 


4 


4.6e-08 


1.2e-08 


3.3 


2 


9e-09 


4.0 


3 


6e-08 


14 


l.le-08 


3 


4 


3.8e-07 


l.le-08 


3.3 


3 


7e-08 


4.0 


5 


6e-07 


15 


7.3e-08 


3 


4 


2.5e-06 


7.3e-08 


3.3 


2 


3e-07 


3.8 


1 


Oe-06 


16 


8.8e-08 


3 


4 


8.8e-06 


8.8C-08 


3.3 


5 


2c-07 


4.0 


1 


4c-06 


17 


7.6e-06 


3 


5 


4.5e-05 


7.6C-06 


3.4 


4 


le-06 


4.0 


3 


2c-05 


18 


2.2e-05 


3 


4 


2.7e-04 


2.2C-05 


3.3 


2 


6c-05 


3.8 


2 


4c-04 


19 


1.2e-04 


3 


4 


6.1e-04 


1.2e-04 


3.3 


1 


lc-04 


3.8 


2 


8e-04 


20 


9.4e-04 


3 


4 


4.7e-03 


9.4e-04 


3.4 


4 


7e-04 


3.9 


4 


5c-03 



on the accuracy of roots computed by qs_dqds. Structured QR iterations algorithm 
compan_qr was the least accurate on this test. 

Table 7.2 
Reversed Wilkinson's first polynomial. 



with balancing without balancing 







qs_dqd 


s 


roots 




qs_dqd 


s 




compan_qr 


roots 


n 




£ 


ni 




£ 




p 


ni 




ni 


£ 


10 


1 


6e-10 


3.3 


2 


4e-09 


1 


6e-10 


3 


3 


3.0e-07 


4.8 


6.7e-07 


11 


1 


le-10 


3.3 


1 


8e-08 


1 


le-10 


3 


3 


1.5e-04 


4.2 


1.5e-04 


12 


1 


6e-09 


3.3 


1 


8e-07 


1 


6e-09 


3 


3 


2.7e-02 


4.8 


1.3e-04 


13 


3 


Oe-09 


3.2 


1 


4e-06 


3 


Oe-09 


3 


2 


l.le-01 


7.8 


1.4e-01 


14 


5 


5e-08 


3.3 


1 


Oe-05 


5 


5e-08 


3 


3 


1.3e-01 


5.7 


2.4e-01 


15 


6 


6e-08 


3.3 


1 


9e-04 


6 


6e-08 


3 


3 


4.6e-01 


7.0 


5.1e-01 


16 


2 


Oe-06 


3.2 


1 


Oe-03 


2 


Oe-06 


3 


3 


4.8e-01 


6.9 


6.8e-01 


17 


6 


5e-06 


3.2 


1 


3e-02 


6 


5e-06 


3 


2 


l.le-t-00 


7.8 


7.3e-01 


18 


5 


3c-05 


3.2 


1 


6e-01 


5 


3e-05 


3 


3 


l.le-t-00 


7.1 


l.le-fOO 


19 


1 


5e-04 


3.2 


1 


8e-01 


1 


5e-04 


3 


2 


1.4e+00 


10.3 


1.4e+00 


20 


3 


7c-03 


3.2 


2 


7e-01 


3 


7e-03 


3 


3 


1.3e+00 


7.2 


1.9e-h00 



Test 3: Wilkinson's second polynomial p{x) = 11 (-^ ~ 0.6*). 

1=1 

In contrast to the Wilkinson's first polynomial case, roots of this polynomial are stable 
with respect to small perturbations of coefficients, although the ratio of the largest to 
the smallest root grows exponentially. Results of numerical tests with this polynomial 
are presented in Table 17731 Both roots and compan_qr algorithms failed to compute 
roots accurately for polynomials of degrees higher than 30. In the case of roots even 
balancing did not improve the situation. On the contrary qs_dqds maintained high 
accuracy of the computed roots for polynomials of degrees up to 50. 
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Table 7.3 
Wilkinson's second polynomial. 



with balancing without balancing 





qs_dqd 


s 


roots 


qs_dqd 


s 


compan_qr 


roots 


n 


£ 


ni 


£ 


£ 


ni 


£ 




ni 


£ 


10 


4.6e-13 


2.7 


3.7e-13 


4.8e-14 


2.8 


1.8e 


-03 


3.3 


1.7e-03 


20 


8.2e-13 


2.4 


2.8e-09 


6.4e-14 


2.5 


2.1e^ 


-00 


7.8 


3.9e+00 


30 


8.1e-13 


2.0 


3.2e-04 


2.1e-13 


2.2 


1.3e^ 


-04 


11.4 


6.9e+05 


40 


8.4e-13 


1.8 


2.4e+00 


1.8C-13 


2.0 


2.1e^ 


-07 


7.7 


1.7e-h08 


50 


1.6e-ll 


1.6 


1.3e+01 


2.5e-13 


1.9 


7.6e^ 


-09 


7.6 


3.4e-l-10 



Test 4: To demonstrate the robustness of the proposed algorithm we conducted 
1000 random tests with polynomials having roots randomly distributed as X ■ 10^^, 
where X and Y are random variables uniformly distributed on [—1,1]. Such polyno- 
mials have both positive and negative roots of very different magnitudes. We report 
mean and standard deviation of the relative error in the tests in Table Ol Tests 
indicate that qs_dqds becomes more accurate (on average) than roots algorithm as 
polynomial's degree increases. 

Table 7.4 

Polynomial with roots distributed as X ■ 10^^, X, "K ~ f/[— 1, 1] . 



with balancing without balancing 







qs-dqds 


roots 


qs_dqds 


roots 


n 




£ ni 


£ 


£ ni 


£ 


10 


e- 


13±e-12 2.2±0.3 


e-12±e-ll 


e-09±e-08 2.3±0.3 


e-K07±e-l-08 


20 


e- 


12±e-ll 2.4±0.2 


e-09±e-08 


e-08±e-07 2.4±0.2 


e-l-07±e-h08 


30 


e- 


ll±e-10 2.4±0.2 


e-07±e-06 


e-08±e-07 2.4±0.2 


e+06±e+07 



Test 5: Algorithm qs_dqds can be applied to polynomials represented in a basis of 
arbitrary orthogonal polynomials. We confirm this by presenting results of numerical 
tests with Wilkinson's polynomial represented in the basis of Chebyshev polynomials 
of the second kind. Chebyshev polynomials were scaled to the interval [a — 1,6-1-1], 
where a and b is the smallest and the largest root of the polynomial of interest, 
respectively. We compared performance of the new algorithm with MATLAB's eig 
routine. Results are presented in Table FTSl 

Conclusion. We extended dqds algorithm of Fernando and Parlett [16l [20] from 
tridiagonal to quasiseparable matrices only in the Hessenberg case but qds algorithm 
to all quasiseparable matrices. New algorithms use just 0{n) operations per iteration. 
Hessenberg quasiseparable matrices include, among others, important companion and 
comrade matrices and, hence, can be used to find roots of polynomials represented in 
different polynomial bases, e.g. monomials, orthogonal polynomials. Our findings are 
supported by preliminary numerical experiments with polynomials having real roots. 
In our implementation we used naive shift strategy and single shifts only. There 
are numerous extensions and improvements of the algorithm left for future research 
including the development of more sophisticated shift strategies for better robustness 
and stability as well as double shifts for handling complex roots. 
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Table 7.5 

Wilkinson's first polynomial in Chebyshev basis. 



with balancing without balancing 





qs_dqd 


s 




■ 


qs_dqd 


s 




■ 


n 


s 


ni 


e 




£ 


ni 


£ 




10 


6.5e-13 


4.1 


4.4e- 


15 


6.5e-13 


4.1 


1.9e- 


15 


12 


1.5e-12 


3.9 


6.5e- 


15 


1.5e-12 


3.9 


3.4e- 


15 


14 


5.9e-13 


3.9 


1.6e- 


13 


5.9e-13 
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